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Thermal field theory is applied to particle production rates in inflationary models, leading to 
new results for catalysed or two-stage decay, where massive fields act as decay channels for the 
production of light fields. A numerical investigation of the Boltzmann equation in an expanding 
universe shows that the particle distributions produced during small amplitude inflaton oscillations 
or even alongside slowly moving inflaton flelds can thermalise. 
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^ , ' Inflationary models give a picture of the early universe that has shown spectacular agreement with observation 
O! [11, Si- AH inflationary models require a mechanism for reheating the universe to take it from the vacuum dominated 
inflationary phase to the hot radiation-dominated universe which we know must follow. The details of particular 
reheating mechanisms depend on the interactions between the inflaton and other fields, but the underlying process is 
particle production which fills the universe with radiation. 

The original work on reheating in the 1980's introduced particle production in an ad hoc fashion, assuming the rate 
1—1 of particle production by the inflaton cj) in the limit of small 4' was proportional to 0" 0, H 0] ■ Most authors settled 
for n = 2, which has the advantage of being equivalent to a simple friction term in the inflaton equation of motion. At 
^p^, around the same time, a less ad hoc approach was based on particle production caused by an inflaton field oscillating 
^ i' about the minimum of the inflaton potential after the end of inflation 0, . This latter approach appeared to be the 
D ] more consistent, and it is still widely used today. 

^ • Following this early work, there where several attempts to apply new ideas in thermal field theory to the reheating 
[ process. These usually focused on finding effective field equations for the infiaton. The small-0 equation of motion 
PsJ ■ was derived first, using linear response theory 0, This has been used in the theory of warm inflation (Til. [T^. [T3| . 
^ ' More general forms of the inflaton field equation, which where not limited to small time derivatives and could be 
On . applied to the oscillating inflaton, followed later [3, [lB| ■ 

' Renewed interest in particle production was i gnit ed b y th e discovery of preheating, a nonperturbative process of 



" inflaton decay though parametric resonance l3, Tsl [T9I Hol [2l|. Like the earlier work, preheating involves particle 
production from an oscillating inflaton field. Preheating is a result of very large amplitude oscillations in the inflaton 
field. Large amplitude oscillations are a feature, though not necessarily a desirable one, of most single-field inflationary 
models. 



Preheating may or may not occur depending on the details of the particle model. In this paper we shall focus 



00 

• • , on models with a mechanism which we call catalysed or two-stage decay 22|. Unlike in the case of pre-heating, we 
. ^ " examine what happens when the fields coupled to the inflaton are too massive to be produced directly. Instead, the 
k> , fields can act as decay channels for the production of light fields. This requires the kind of mass hierarchy which is 

; I possible in supersymmetric models. Examples include Grand Unified extensions of the standard model. 

Ci ^ In section II we obtain formulae for the particle production rate due to changing particle masses or couplings 
and show that to a reasonable degree of accuracy the expansion of the universe can be neglected in making these 
calculations. In section III we consider several models of particle production that could arise during infiation, starting 
with the oscillating and slowly-evolving models, for which particle production rates are well-known and easily checked, 
and then we apply our method to the catalysed decay. 

In section IV we use the production rates calculated in the previous section to take a closer look at the evolution 
of a system which may be expanding and departing from thermal equilibrium. Thermalisation has previously been 
considered in the context of reheating through numerical solution of classical non-linear field theory psi . [2^ or using 
a numerical solution for the non-equilibrium particle propagators psi . [26j . We solve the Boltzmann equation in an 
expanding universe with a source term representing particle production. This allows us to consider models which 
would be too difficult to analyse directly using numerical approaches to quantum field theory. 
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II. PARTICLE PRODUCTION 



In an ideal situation we would like to track the formation of particles during and after the inflationary era to give a 
full description of the reheating process. This requires a workable definition of particle number. This definition need 
not be unique, but at least it should agree with the usual definition of particle number after inflation has ended. We 
shall use a definition of particle number which introduces a free field which instantaneously has the same field value 
and momentum as the interacting particle field. 

A fundamental problem we face is that the particle production might not have a local description. However, we 
might hope that simple situations occur where the particles are produced at a rate depending only on local conditions, 
for example local field values or temperatures. For this reason, we focus on particle production rates. 

We shall consider the production of particles due to a background time-dependent inflaton field (f>. The particles 
will be associated with a field a. 



Particle number 



Following Morikawa and Sasaki [lO|, we define particle creation and anihilation operators for a fiducial free field 
which coincides with the interacting field a and momentum tt at time 



a^{p,t) = ujpd-{-p,t) - iTr{p,t) 
a{p, t) — Ljpd-{p, t) + in{-p, t) 



where ujp = {p^ + m^)^/^ may depend on time, and 



a{p,t) — J d^xa{x,t)e 



-ip-x 



(1) 

(2) 



(3) 



The local number density is assumed to be spatially homogeneous. The number density in phase space n[p) can then 
be defined in terms of an ensemble average using the fiducial free field. 



1 

2uj, 



■{a\pi,t)a{p2,t)) - {27rf S{pi ~P2)n{pi). 



We shall express the density function in terms of the Wightman function G2i{p, ti, 12), defined by 

(0-(pi,ti)(T(p2,<2)> = (27r)3(5(pi -p2)G2l(pi,tl,t2)- 

The expression for the density function becomes 



n{p, t) = 



1 



- — {ujp - idti){ujp + idt^)G2i{p,ti,t2) 



(4) 



(5) 



(6) 



where [. . .] will be used to indicate when a function is to be evaluated at t\ =t2 = t. The frequency Up always refers 
to the value at time t, unless we state otherwise. 



B. Particle creation 



In order to use Eq. ([6]), we would have to solve field equations for the Wightman function with suitable initial data, 
giving a particle density which is typically a non-local function depending on the history of the background field. 
Instead of working directly with the density function directly, we shall look for a local approximation to the particle 
production rate. 

After some elementary manipulation, the particle production rate in phase space obtained by taking the time 
derivative of Eq. ^ becomes 



fl — ^mass ^int . 

The first term represents particle production due to the changing particle mass, 

1 



2uj„ 



(ujp - idtj + (ujp + idt^) (up- idt^){ujp -f idt^) G'2i(p, ti, t2) 



(7) 



(8) 
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The second term represents particle production at fixed particle mass, 

(9) 



^iiit 



p 



We can put these expressions into a more useful form by introducing the self-energy. Since we are interested in 
the evolution of operators in a given initial state it proves convenient to use an 'in-in' formalism, and we use the 
Schwinger-Keldysh version (27l . [28j . Propagators carry two extra internal indices a and 6, where the indices a and b 
take the values 1 or 2. Following Calzetta and Hu [23|, we raise and lower indices with a metric Cab = diag(+l, —1). 

The Schwinger-Dyson equations for the in-in formalism read [2^ . 



{dl+L0l{t2))G2l{pMM) = - J dt'G2''{p,h,t')Eal{p,t',t2) (10) 
{dl+iul{h))G2lip,h,t2) = - f dt'^2''ip,ti,t')Galip,t',t2). (11) 



For the terms in the particle production rate which have just one time derivative of the propagator, we use the LSZ 
trick of introducing a time integral 



(a*, -zc.p)G2i(p,ii,i2)|,,=, = / dt2e-^^''('~'-\dl+cul)G2i{p,h,t2) (12) 
{dt,+tiUp)G2i{p,h,t2)\,^^, = / dhe-''^''^'^-'\dl+cul)G2i{p,h,t2). (13) 



We can now express the particle production rate in terms of integrals of the propagator and the self-energy which 
have a suitable form for applying perturbation theory. 

We shall consider the leading order in perturbation theory. First of all, let 

ujl=p^ + ml + c/(^Ht) (14) 

where is a constant and (j){t) is given. Suppose also that S = 0{g^). This is the kind of situation would arise, for 
example, given an inflaton and an interaction Lagrangian density C — g^cji^a"^ / A:. 
The leading order result, using Eq. ^ and Eqs. (jlOlllgp . is that 

/ dt2-- {^\t)~4>^{t2))G2l{p,t,t2)\ (15) 

J-oo 2Wp J 

where G2i{p, t, <2) is a free Wightman function for the ct- field with the shifted mass (which need not be in the vacuum 
state). 

For the next term, we require the derivative of the Schwinger-Dyson equation JTI 



{dl +Lul){dl +Lul)G2l{p,h,t2)^l^2l{p,tut2) + Oig'^). (16) 

Use this together with Eqs. jS]), ^ and (fTS]) . 

nint = Im|2y dt2 — E2i(p,t,i2)j . (17) 

where S21 is the self-energy of the cr-field at leading order. We can see from this expression that riint is the part of 
the production rate which is associated with the imaginary part of the self-energy. 



C. Curved space 



The formulae for the particle production rates found in the previous section neglected the expansion of the universe. 
In this section we shall seek to show that expansion can be neglected to a reasonable degree of accuracy when particle 
momenta are larger than the expansion rate. 
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Consider a spatially flat universe with scale factor a and constant expansion rate H. Solutions to the wave equation 
in de Sitter space can be decomposed into particle modes 30] with comoving wave number k which satsify 

{d^ + 3Hdt+ujl)f{k,t) = (18) 

where 

= k'^/a^ + + 12^H'^. (19) 
A suitable normalisation is to use //* — /*/ = i/a^. The modes can be expressed in terms of Hankel functions, 

f{k,t)^^Hk-'/'.y'Hl'Hz), (20) 

where z = k/{Ha) and v"^ = 9/4 - rr? jH"^ - 12^. 

Consider a non-interacting field with Wightman function 

G2i(fc, txM) = (iV(fc) + 1) /(fc, ii)/*(fc, t2) + iV(fc)/(fc, t2)/*(A, tl). (21) 

The phase space number density is defined as before, 



n(p, t) 



(22) 



This is a function of the momentum p = k/a of a locally defined flat-spacetime theory. The factor c? is present because 
the Wightman function has undergone a shift in normalisation due to the use of comoving modes. We substitute the 
Wightman function (PTjl and use the modes (|^D|) . The result is quite complicated in general, but the main features 
can be seen in the conformal case m = and ^ = 1/6. In this case, 

n = UdeS + "rod- (23) 

where 

Ifc / 3ffa\ 1 
n,.soJ, = --(^l + -_j--a;, (24) 

k( 3Ha\ 



a 



The non-vanishing contribution to the density function in the de Sitter vacuum is a reflection of the fact that the 
particle number was defined using the physical momentum. The particle number defined this way is analogous to the 
response of a particle detector. We can see from Eqs. (|24l) and (|25|) that n Ri N for p >> H, so that in this limit we 
recover the flat space results. 

Similar considerations apply also to the particle production rates. However, it is important to bear in mind when 
calculating particle production rates that -h is evaluated at constant p and not constant fc, 

^ -HpUr-] + Ur ■ (26) 



9t J p \dp J t V 

The first term on the right of this equation represents the reduction in particle density caused by the expansion of 
the universe. The second term, representing the particle production, goes over to the fiat space-time results when 
p^ H. The redshift term is analysed further in Sect. IIVI 

It is interesting to integrate Eqs. p4|) and (|25|) to get a formula for the energy density. 



deS 



(27r)3 a ^ ' \ 2 k 



The vacuum energy density of de Sitter pdes space comes from integrating the left hand side of Eq. ((24)) . after we 
apply suitable regularisation methods. The full expression combines both curved space and thermal effects, with 
thermal effects dominating the integral for p >> H. 



5 



III. EXAMPLES 



We shall take a closer look at four different models and calculate some particle production rates using the formalism 
described in the previous section. These models are typical of what one might expect in the context of inflation. Some 
of these results have been obtained before using other methods, and these are included to check the consistency of 
the new approach. 



A. Oscillating fields 

The first example we consider is the particle production from small amplitude oscillations of an inflaton or other 
background field. In many inflationary models, this type of particle production would be eclipsed by preheating from 
large amplitude oscillations. However, this is not always the case, so that even this simple example may be of interest. 

The background field we take has 

= 0o(l + ecos TO^i). (28) 

where (j)o is the stable vacuum value of the field and e is small variable which varies slowly on the oscillation timescale. 
The field (f> is coupled to a field cr with effective frequency LUp, 

ul^p' + ml+g^4>^-g^4>l. (29) 

In this case we have introduced a shift so that the mass is mo- when (jj = (/jq. 

This particle production problem was solved long ago (31i] . The total 2— particle production rate is given by the 
standard formula for particle decay, 

"■^^ (30) 

where the momentum p is determined by momentum conservation and the reduced matrix element M is defined by 

(Pi, P2IO) = X(27r)*(5(pi + P2)(5(cjp, + ~ m^) (31) 
To leading order, the interaction with the background through Eq. ([^5]) gives 

(\ 1/2 
1 f > 2m,. (32) 

J 

The dependence on the scalar field energy density = 0Qe^m^ is usually factored to define the reheating coefficient 

rby, 

r = ^«5VL_ (33) 

Since g is typically very small, this type of perturbative reheating is quite inefficient and would take several Hubble 
times to complete. 

We can consider the same problem, but using the general result for the time derivative of the density function ([H]). 
The free Wightman function for the vacuum state is given by 

G2i{t) = -^e-"^"'. (34) 



2ujp 



After integration over time, we have 



^ = 75 ^0^ -f sin i^^i) ^i^P - "^0/2) + Tff ^0^ -J ,2_^2v (35) 
The total particle creation rate integrated over momentum is 

1-^1 sin2(m^t) + ^g^l^ "^^' sin(2m^t). (36) 

This result appears complicated, but this is due to the presence of transient terms. Such terms are to be expected when 
we try to calculate the particle production instantaneously. Over several oscillatory cycles, however, the production 
rate averages out and we recover the scattering theory result p2p . The particle production in this case is effectively 
localised as long as we consider times longer than the oscillatory cycles. 
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B. Derivative expansions 



Another situation where we can have locahsed particle production rates is in the 'adiabatic' hmit, when the inflaton 
has a smah time derivative. This is a speciahsed form of particle production which does not usually occur at leading 
order in perturbation theory. An important exception occurs when the system starts out and remains close to thermal 
equilibrium. This type of particle production was first discovered by Hosoya and Sakagama @ and by Morikawa and 
Sasaki 

We start again from the general result for the particle creation rate ([8|), using the adiabatic approximation 



'-2(^2) = <f>-'{t) ~ ^\t2) ~ 20(O0W(i2 - t). (37) 



We introduce transforms, 



After integration, we arrive at 



G2iit-h) = I ^e-'^-(*-*^)G2iH (38) 

- / ^e-(*-*^)5^2(c.) (39) 
, Zn 



'■9k^Z^, (40) 



It only remains to give a formula for the thermal Wightman function. This can be expressed in terms of a spectral 
function p and the thermal distribution function n (for example, see [32j). 

G21 = -i(l + n)p. (41) 

The spectral function typically has a Breit- Wigner form [ssj , 

{cj^ -uj^- 2iujT-^)-^ - (cj2 + 2iujT-^)-^ (42) 

where r is known as the relaxation time. Inserting the Wightman function into Eq. (|40p gives 



h = -g^c^^^^ 43 

This formula can also be derived using the methods of Ref. , and it is closely related to the work of Ref. @ . 

Note that the particle production is exponentially small for temperatures less than the cr— particle mass. In fact, 
Eq. (j40[) vanishes at zero temperature due to a general property of the Wightman function. The best way to view 
this type of particle production is as a type of transport phenomenon, similar to thermal or electrical conductivity. 
Particles are produced as the system responds to the disturbance of thermal equilibrium caused by the changing mass. 
Increasing the relaxation time r gives more time for the mass to change, driving the system further from equilibrium 
and increasing the particle production. 

The particle production takes energy from the inflaton fleld, and we can ensure energy balance by introducing a 
friction term T(f> into the inflaton equation. The total radiation energy is 

riLUp (44) 



(27r) 



The time variation of pr contains a term from the time variation of the ujp, which relates to the change with time of 
the inflaton effective potential, and a term, where 



(27r)3 



,2- 10^.-1^ ■ 



The friction coefflcient obtained from the particle production formula agrees with the friction coefficient obtained in 
the inflaton equation of motion using linear response theory Q. 



7 



C. Catalysed or two-stage decay 



In the third model of particle creation an oscillating inflaton decays into light thermal scalar particles through 
an intermediate virtual boson. This is a natural set-up, because many particle models contain both heavy and 
light fields, with the light particle masses protected by supersymmetry. Particles which couple to the inflaton will 
tend to be massive, and may well be too heavy to be produced directly by the inflaton oscillations. Preheating is 
supressed, and we have to rely on perturbative particle production effects. We shall suppose that the model is part 
of a supersymmetric theory, which protects the flatness of the inflaton potential and the masses of the light fields, 
doing away with the need to fine-tune their coupling constants. 

The mass of the light particles in this example is independent of time, and we use the second formula ([9]) for the 
particle production rate. Having a fixed mass removes some of the ambiguities in the definition of the particle number, 
and leads to a 'cleaner' result. 

The heavy field is denoted by x ^-^d the light field by <j. The interaction Lagrangian we shall take is 

Ci = -\g^^^X^-\hma\^^\a\ (46) 

where g, h, X and m are constants. We can choose m = gcpo by redefining h if necessary. The self-interaction term of 
the light fields is included to allow them to come to thermal equilibrium. We shall consider particle production into 
the vacuum and also in the presence of thermal radiation. 



p — f hm hm» — 

FIG. 1: The Feynman diagram contributing to the imaginary part of the a self energy E12. 



The background is as before Eq. (j28p . but now > m^. The first non-trivial contribution to the imaginary part 
self-energy is given by the Feynman diagram show in figure [1] We define the fourrier transform as in the previous 
section and then the diagram contributes 

„ , s A, o o f d?k duji dujy duj-i ,-,/-** ^ „ 

E.rb,M.) = gWj j^^^^^—''~''-'G.^AP-Kt~t,) 

Gx2"(k, c^i)02(^^ _ i03)G^a\K i03)^\^3 " u;2)G^bi (k, c^i), (47) 

where a subscript has been used to distinguish between x and a propagators. 

We concentrate on the low energy spectrum p « m^, when we can use a low momentum approximation for the x 
propagators, 

Gx2'(k,t^) « (48) 



G^2^(k,c.) « ^e{uo) (49) 



G^n(k,u;) « --^ (50) 

where 9{oj) is the heaviside function. The middle equation follows from Eqs. (|¥T|) and ([1^ . where is now the heavy 
particle decay width and a — ALOp/r oc h'^m?. We use the free Wightman fmiction for the a field with occupation 
number n. 
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With these approximations, using Eq. ([9]) for the particle production rate, 



27r2 



-np), 



where 



F{p) 



{mil, - LUp- LUk) (1 + . 



For small 771^, the integral gives a dilogarithm function, 

^(P) = —dilog [e 



{m,-ujj,)/T 



m 



(51) 



(52) 



(53) 



The function F{p) has been plotted in figure O As might be expected, the vacuum particle production rate peaks 
when the momentum equals half the inflaton mass. The physical process behind the particle production this time is 
a decay (f) — > 4(t, using two intermediate virtual x bosons. The four particle decay is reflected in the broad width of 
the peak, compared to the resonance in the 2— particle decay in Sect. IIII Al 



0.2 



0.15 



^„ 0.1 



0.05 



T=0 



0.2 



0.4 0.6 



FIG. 2: The momentum dependence of the particle production rate for the two-stage decay with an oscillating inflaton (model 
C). The function p^F{p) has been plotted for nia = at T = and T — m^. 



In the zero temperature limit, the reheating coefficient F which describes the rate of production of radiation becomes 

P0 480^4 ^8 ■ y ' 

For a ^ h'^m?^ this is smaller than the corresponding result in Sect. IIII Al for the to^ > regime. However, in 
supersymmetric models, the couplings do not have to be especially small. Furthermore, if there are many species of 
light scalar fields, then the the final result scales with the number of fields. 



D. Derivative expansion for catalysed or two-stage decay 



The final example is another 'adiabatic' process, but this time the low temperature behaviour is suppressed by a 
power law instead of the exponential suppression found in Sect. IIIIBI The inflaton decays into light thermal scalar 
particles through an intermediate virtual boson as in the previous example. This model was introduced in the context 
of warm inflation ^22i] . but the set-up can occur quite easily in models which contain both heavy and light particles. 
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The interaction Lagrangian is the one used in the previous section, Eq. (|46|) . We take an initial state to be one of 
thermal radiation with temperature T « m^. How the system might come to thermal equilibrium will be addressed 
in the next section. 

The first non-trivial contribution to the self-energy is again given by the Feynman diagram show in figure [1] and the 
expression (|47p . We can use the condition T « to justify a low momentum approximation for the x propagators 
again, now with 

(55) 

) (56) 

(57) 



and energy, which has been calculated explicitly for non-zero temperatures in Ref. 34 1 . 
We use an adiabatic approximation for the inflaton fields as in Sect. IIII B[ 

5<tP{u) = 2i(j)^{2T:5'{uj)) (58) 

where the primes denote derivatives with respect to uj. With these approximations, the self-energy becomes 

5]2i(p,i,t2)=4//i2!^02 f ^|^^e-Mt-t2)G^^^(p_k,t_i2)(a(c.)[H-n(c.)])". (59) 
Now we can apply formula ([9]) for the particle production rate, using the free thermal propagator for the a field, 













i 

ml 


(SB and El 


with a 



= _4g4/^2^^2 / ^ {n"{-0Jk - ujp){l + n{uk))} . (60) 



^pm\ J (27r)3 LUk 

So far, we have not had to assume a particular form for the distribution function n. However, if n is the thermal 
distribution function, then the integral can be done approximately in the small m„ mass limit. 



flint = -TT^ 8"^(P) (^1) 



where m = g(f> and 

F{p) = |n(c.p) -1 (l - ^/^) | . (62) 

The function F{p) has been plotted in figured Note that the pre- factor is identical with the pre-factor in Eq. ([5T|) for 
the oscillating case if we use the period averaged value of (p^. The momentum distribution is quite different, however, 
and the particle production vanishes as T — > 0. 



The reheating coefiicient F can be obtained by integrating over momentum, using a — AtOp/r « const, 

? 1 2 A 

This agrees with the friction coefficient calculated from the effective field equations in ref. [33| • 

IV. THERMALISATION 

The adiabatic particle production rates where calculated under the assumption that the system was close the thermal 
equilibrium. In this section we shall examine the validity of this assumption by solving the Boltzmann equation in an 
expanding universe. This will also give us an opportunity to consider systems which depart from equilibrium. 
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FIG. 3: The momentum dependence of the particle production rate for smaU (f) with the two-stage decay (model D). The 
function p^F{j)) is plotted. The thermal distribution p^n is shown for comparison. 



We shall adopt a pseudo-particle approximation where the Wightman function take a thermal form, but with an 
arbitrary distribution function n = n{p, t). The particle number will evolve according to 



(64) 



where Sp represents particle production, Sr represents particle dilution due to the expansion of the universe and Sc 
is the Boltzmann collision term 

The particle production rates calculated in Sects. IIII BlUIlDj are still valid with the new distribution functions and 
can be used for Sp. The expansion effect represents a stretching of the physical wavelengths of the modes by the scale 
factor a. This term was evaluated in Eq 



Sr — HpdpTi. (65) 
The collision term for 2^2 particle scattering from the quadratic term in the Lagrangian density Eq. (|46p is 



[ tfp2 d^P3 (Ppi /^_,_5 



2cJo ./ 2cjp2 2cjp3 2cjp^ 



{2ti)-H{P + P2-P3- PA)Bip,pi,p2,P3), 



where P = (ujp, p) and 



B{p,Pi,P2,P3) = (1 + "(p))(l + n{p2))n{p3)n{pi) - n{p)n{p2){l + n{p3)){l + n{p4)). 

This term preserves the total particle number as well as the total energy. 

Multiplying the Boltzmann equation by LOp and integrating gives the total energy equation 



where the source term is 



S = 



d'p ^ 



(66) 



(67) 



(68) 



(69) 



In the oscillating inflaton case S — Tp,j, and in the slowly evolving inflaton limit S — F^^, where expression have been 
given for F in Sects. IIII Bl HITdI Both types of reheating coefRcient can be combined into 



AHpr = T{p^+p^), 



(70) 



where is the averaged pressure term. Various forms of this equation have been used in the past to study reheating 
[d| and warm inflation 13511. 



11 



The simplest way to analyse the Boltzmann equation (|64p is to take a close-to-equilibrium approximation, intro- 
ducing a thermal distribution function nT and defining the effective temperature by matching the energy density with 
the actual distribution function, 

y (2^(^-"tH = 0. (71) 

We can use the thermal particle production rates calculated earlier and introduce a thermal relaxation-time r,. to 
simplify the collision term. The Boltzmann equation we have to solve is then 

RF{p) + Hpdpn-T-\n~nT), (72) 

where R is the prefactor in the particle production rates Eq. ((5T|) or Eq. ((6T|) . We have integrated this equation 
numerically using a fourth order Runge-Kutta scheme for the time derivatives and second order differences for the 
momentum derivatives. This numerical procedure is very fast and stable, with most of the results given below taking 
less than one second on a IGHz processor. 



A. Oscillating phase 

We consider a period of reheating with an oscillating inflaton and >> pr- This regime ends, according to Eq. 
(j70p . when H T. During this period, the pressure averages to zero over the oscillation period and the universe 
expands like a pressure free cosmological model, with 

H{t) = H{0) (^1 + lHiO)?j , (73) 

R{t) = R{0) (^1 + lHiO)?j . (74) 

The second equation follows from R on p^. 

Some numerical results for the momentum distribution obtained from Eq. (|72p are shown in figure |31 The distri- 
bution thermalises, and does so more quickly with smaller relaxation times as might be expected. The momentum 
distribution of the source term shows up clearly at early times, before the relaxation has taken effect. 




FIG. 4: The stationary momentum distribution for the oscillatory phase in the two-stage decay (model D) using the relaxation- 
time approximation. The relaxation times are Tr = 0.1/H{0) (left) and Tr = 0.05/H{0) (right). As might be expected, shorter 
relaxation times produce a spectrum which is closer to thermal equilibrium. The constant -R(O) = 50-ff(0) and = 2H{0). 
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FIG. 5: The time evolution of the effective temperature for the oscillatory phase with the two-stage decay (model D) using the 
relaxation-time approximation. The de Sitter temperature H/2ti is shown for comparison. The relaxation time r,. = 0.1/7^(0), 
the constant _R(0) = 50J/"(0) and = 2ff(0). 



The initial temperature for the numerical solutions has been set equal to the de Sitter temperature H(0)/2Tr, to 
be consistent with the assumptions used in the particle production calculations. The evolution of the temperature is 
shown in figure O After a sharp rise to a maximumjthe temperature falls off as t^^^'^. This agrees very well with the 
analytic solution to the total energy equation (|7(I)) d, [s^ . 



B. Slow-roll phase 

Small values of (p are characteristic of the slow-roll phase of inflation. The particle production rates calculated 
in sect IIII Dl can be applied to the slow-roll phase, provided we can justify the thermal hypothesis which was used. 
During the slow- roll phase of inflation, both H and F vary very little over several Hubble times, and we can treat 
them as constants in the Boltzmann equation (|72p . 




FIG. 6: The stationary momentum distribution for different relaxation times in the two-stage decay (model D) using the 
relaxation-time approximation. The relaxation times are — 0.2/H (left) and Tr — 0.05/H (right). As might be expected, 
shorter relaxation times produce a spectrum which is closer to thermal equilibrium. The constant R = 15H and nia ~ 0.25H. 
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FIG. 7: The time evolution of the effective temperature for different initial conditions with the two-stage decay (model D) 
using the relaxation-time approximation. In each case, the momentum distribution reaches a stationary state with the effective 
temperature shown. The relaxation time Tr = 0.1/ H , the constant R = 15H and ma- — 0.25H. 

Numerical solutions for two different parameter sets are shown in Fig. [6l These show the existence of an attractor 
with non-zero temperature and a spectrum close to thermal equilibrium. The final momentum distribution does not 
show any dependence on the initial distribution, but it is dependent on the relaxation time, a shown in figure[5] Small 
relaxation times, corresponding to relatively large values of the self-coupling A, lead to nearly thermal spectra. 

The parameters for the numerical solution where chosen to place the temperature in the range T > H required for 
consistency of the particle production calculations. The time evolution of the temperature shown in figure [7] agrees 
very well with the analytic solution to Eq. ([70]) when T cx T'^ (see Eq. which has the form 

T = To, (1 - e-^*) . (75) 

This shows clearly how the expansion of the inflationary universe need not lead to a supercooled state when particle 
production is taken into account. 

C. Thermalisation with the full boltzmann collision integral 

In the above work we have introduced the thermal relaxation-time to approximate the thermalisation effects of the 
Boltzmann collision term. We can check the validity of this approximation by solving the the Bolzman equation with 
the full 2^2 particle scattering term ((66)) . Following the work of Refs. [siI-IstI. we can eliminate the delta-functions 
and reduce the integral from 9 to 2 dimensions, which gives 

D f 

5c = / 9{ujp, - nic)nim{p,p2,P3,P4)B{p,p2,P3,P4)dujp^dcjp^, (76) 

WpP J 

where D — A^/647r^ and uj^^ = (p, , ) is obtained from energy conservation. 

We have solved equation (|7^ numerically with the new expression for Sc, focussing on the two stage decay model 
(model D). Again, we used a fourth order Runge-Kutta scheme for time derivatives and second order differences 
for the momentum derivatives. The collision term was evaluated using a 2D Simpson's rule integrator. In order 
to remove instability problems at low momenta we damped the source term with a factor p^Hjp^ + H'^), which is 
consistent with our calculation of the source term which cannot be used for p less than H. We also avoided using a 
very fine momentum mesh that would bring in grid points at very low momentum. Solving with the full collision term 
is computationally far more demanding than using a relaxation-time approximation. For reasonable mesh sizes the 
total integration times are approximately an hour on GHz processors, compared to one second for the relaxation-time 
approximation. 
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Numerical results for the full collision term with two stage decay model are shown in figures [8] and [9l obtained 
using the same values for constants R and as before. The distribution reaches a stable non-zero temperature as 
expected and is consistent with the findings using the relaxation-time approach. 

Comparison of Fig. [S] and Fig. [5] suggests that the relaxation-time = Q.IH^^ corresponds to 13 « 10. This 
example is strongly self-coupled. However, it is possible to argue that value of D needed for thermalisation decreases if 
we increase the particle production rate. According to dimensional analysis, the relaxation time should be proportional 
to the inverse temperature. The numerical example has T = 2H, hence D « 2/(TTr). We therefore predict a similar 
distribution function to Fig. [5] for D < 1 when the particle production is increased to give an effective temperature 
T > 20H. 
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FIG. 8: The stationary momentum distribution in the two-stage decay (model D) using the full collision integral gives a check 
for consistency of the relaxation time approximation used in Fig. [6] The parameters are R = 15H and rricr — 0.25// and 
D = 10. The plot is comparable to Fig. |6]with a relaxation time Tr = Q.l/H. 
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FIG. 9; The time evolution of the effective temperature for different initial conditions with the two-stage decay (model D) 
using the full collision integral. As with the relaxation time approximation, the momentum distribution reaches a stationary 
state with the effective temperature shown. The constant D = 10, with R = 15// and mo- = 0.25//. 
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V. CONCLUSION 

We have attempted to produce a uniform description of particle production during the early universe which can 
cope with oscillating and slowly varying inflaton background fields. We have concentrated mainly on a two-stage 
decay process where the inflaton decays into light radiation fields through an intermediate heavy boson. 

Thermalisation of the particles has been described by solving the Boltzmann equation in an expanding universe. We 
have found that the thermisation and particle production can be combined to produce a prediction for the momentum 
distribution in the radiation fields. In many cases, where the self-interactions allow, the distribution approaches a 
thermal distribution. 

Our results for particle creation and thermalisation in the case of a slowly-evolving inflaton field are fully consistent 
with the thermal dissipation processes predicted in warm inflation 33, 3i|. Most of these models have have assumed 
that the radiation remains close to thermal equilibrium, and we have found that this occurs when the the self-coulping 
of the radiation field is sufficiently large. 

The particle production rates, and therefore the thermal dissipation rates, are still significant even when the the 
distribution function departs substatially from thermal equilibrium. Distortions to the spectra due to the finite 
relaxation time of the radiation may have observational consequences if the thermal fluctuations are the source of 
density fluctuations in the cosmic microwave background [III Sj US • Further work would be worthwhile to find 
the effect this may have on the spectrum and as a source of non-gaussianity. 

The reason for considering the two-stage decay lay partly in the fact that there are light fields whose masses are 
protected by supersymmetry. In a supersymmetric model, the bosonic decays which we have considered would be 
accompanied by fermionic decays. The extension of the present results to fermions is tedious, but straightforward. 
Fortunately, fermion decays tend to be suppressed at low temperatures, compared to the bosonic ones 34], and so it 
should be reasonable to ignore them. 

An interesting regime occurs when the temperature is comparable to the expansion rate. Both the thermal equi- 
librium and flat-spacetime approximations break down in this limit. We have suggested ways to deal with this case 
using curved space methods in Sect. Ill C[ but further work along these lines would be of interest. 
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